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Abstract 

A first-principles molecular dynamics (MD) scheme is presented on the basis 
of the density-functional (DF) theory with use of the the quantal hyper-netted 
chain (QHNC) approximation. The DF theory brings about exact expressions 
for the ion-electron and ion-ion radial distribution functions (RDF's) of an 
electron-ion mixture as a model of a simple liquid metal. These exact expres- 
sions prove that an ion-electron mixture can be treated as a one-component 
liquid interacting only via a pairwise interaction in the evaluation of the ion- 
ion RDF, and provide a set of integral equations: one is an exact integral 
equation for the ion-ion RDF and another for an effective ion-ion interaction, 
which depends on the ion configuration specified by the ion-ion RDF. Hence, 
after some approximations are introduced, the MD simulation can be per- 
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formed to get the ion-ion RDF using the ion-ion interaction determined so as 
to be consistent to the ion-ion RDF: the MD simulation and the procedure to 
determine the effective interaction from the QHNC equation are performed 
iteratively. This MD simulation coupled with the QHNC equation (QHNC- 
MD method) for the effective interaction provides a first-principles calculation 
of structures of simple liquid metal: the ion-ion and electron-ion RDF's, the 
charge distributions of an ion and a pseudoatom, the effective ion-ion interac- 
tion and the ion-ion bridge function are evaluated in a self-consistent manner 
from the atomic number as the only input. 

We have applied this QHNC-MD method to Li, Na, K, Rb and Cs near the 
melting temperature using upto 16000 particles for the MD simulation. It is 
found that the convergence of the effective ion- ion interaction is fast enough for 
practical application to alkali metals; two MD runs are enough for convergence 
within accuracy of 3 to 4 digits, if the initial effective potential is properly 
set up. The structure factors, thus obtained, show excellent agreements with 

experimental data observed by X-ray and/or neutron scattering. 
61.25.Mv, 61.20.Gy, 61.20.Ja, 71.50.+t, 31.20.Rx 



2 



Typeset using REVTp^X 



I. INTRODUCTION 



The liquid alkali metals have been studied extensively in both experimental and theo- 
retical sides. They can be easily used to test a theoretical approach as the first step, since 
they constitute "simple" metals and "simple" liquids: furthermore, there exist many reli- 
able experimental results to be compared. In the standard theory, a liquid metal is treated 
as a one-component liquid interacting via a binary effective potential, which is determined 
by the pseudopotential formalism; a pseudopotential is introduced either by first-principles 
calculations or by adjusting parameters involved in model potentials to some experimental 
results. In this treatment, the ionic structures are determined independently of the electronic 
structures in a liquid metal. 

It is only recent that a liquid metal is thought of as an electron-ion mixture and the 
ionic structures are determined in a coupled manner with the electronic structures. One of 
such approaches is the Car-Parrinello molecular-dynamics (CP-MD) technique IJ, where a 
liquid metal is taken as a binary mixture of ions and electrons. In the CP-MD method, the 
electron-ion interaction is described by a pseudopotential to produce pseudo-wavefunctions 
which can be accurately represented by small number of plain waves. The CP-MD method 
possesses advantage to avoid the difficult task of constructing an effective ion-ion potential 
required to perform the molecular-dynamics simulation, and afford to provide ab initio 
calculations of metallic systems in principle. However, the most serious problem in this 
approach is that the number of particles used in the simulations cannot be taken large 
a nd a total of time steps performed in the simulations is limited to a small size within 
the present computational resources. In the present work, we propose another scheme of 
first principles molecular-dynamics simulation based on the density-functional (DF) method 
applied to the ion-electron mixture; this simulation method can be performed on the large 
number of particles to last a large size of time steps, since this scheme reduces the electron-i 
on problem to a usual classical MD coupled with a set of integral equations determining an 
effective ion-ion interaction: a problem to determine the ion-configuration structure and the 
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ion-ion interaction in a self-consistent manner. 

Previously, we have proposed a set of integral equations for radial distribution functions 
(RDF's) in an electron-ion mixture on the basis of the DF theory in the quantal hyper-netted 
(QHNC) approximation |2|||. In this QHNC formalism, the bare electron-ion interaction 
v e \(r) and the ionic charge Z\ are determined self-consistently by regarding a liquid metal 
as a mixture of nuclei and electrons (§]. Already, we have applied this approach to liquid 
metallic hydrogen 0, lithium ||, sodium ||, potassium [Q and aluminum || obtaining 
ion-ion structure factors in excellent agreements with experiments. The QHNC equations 
are derived from exact expressions for the electron-ion and ion-ion RDF's in an electron-ion 
mixture: these exact expressions are only formal results derived from the DF theory. A 
molecular-dynamics scheme to treat an ion-electron mixture can be set up on the basis of 
these exact relations, which states that the electron-ion mixture can be regarded as a quasi- 
one component liquid interacting via a pairwise interaction in the description of the ion-ion 
RDF || (hereafter, referred to as the QHNC-MD method). Since the QHNC formalism is 
derived on the electron-ion model where the bound-electrons forming an ion in a liquid metal 
is assumed to be clearly distinguished from the conduction electrons and the overlap of the 
core electrons is negligible, the QHNC-MD method is only applicable to simple liquid metals: 
its application to liquid alkali metals is taken as an ideal test of the QHNC-MD method. 
Thus, the QHNC-MD method has been shown to yield structure factors of liquid alkali 
metals in excellent agreement with experiments as the result of a first-principles calculation 
in the present work. 

In Sec. |H|, we sketch the QHNC formulation: exact expressions for RDF's in an electron- 
ion mixture are obtained from the DF method []10|1 , and the nucleus-electron model is shown 



to provide a bare electron-ion interaction, which should be determined self-consistently. The 
procedure to perform the MD simulation based on the QHNC theory is shown in Sec. [Til]: in 
the QHNC formulation the effective ion-ion interaction used in the MD simulation depends 
on the ionic structure specified by the ion-ion RDF. Therefore, in the application of this 
MD scheme it is important to extrapolate the MD RDF beyond the truncation radius of the 



simulation correctly so as to be used in the determination of effective ion-ion interaction; this 
is exemplified by the method described in Sec. [TTT]. Numerical procedure of the QHNC-MD 
method and the results of its application to alkali liquid metals are described in Sec. [IV]. The 
last section is devoted to a discussion, where the advantage/disadvantage of the QHNC-MD 
method against the CP-MD method is argued also. 



II. QUANTAL HYPER-NETTED CHAIN THEORY 

A simple liquid metal can be thought of as a binary mixture of ions with a definite ionic 
charge Z\ and the conduction electrons; the interactions between particles [i,j =1 or 

e] are taken as pair-wise. The ions constitute a classical fluid, while the conduction electrons 
form a quantum fluid. Let us refer to this mixture as the ion-electron model for a liquid metal. 
Since the ions are regarded as classical particles in the electron-ion model, the ion-ion and 
electron-ion RDF's become identical with the ion- and electron-density distributions around 
a fixed ion in the mixture, respectively ||. Because a fixed ion causes external potentials 
acting on ions and electrons in the homogeneous mixture, the DF theory can give the exact 
expressions for the ion- and electron-density distributions, rii(r \ I) and n e (r\ I), in terms of 
those of noninteracting systems n°(r) under effective external potentials U\ ? {r) [i=l, e] 

Uf(r)=v a (r) + -^-»T t (1) 
oni(r\ I) 

with the use of J 7 -^ and /^ nt , the interaction part of the intrinsic free-energy and the chemical 
potential, respectively 0. As a result, the DF theory provides exact, but formal, expressions 
for the ion-ion RDF g u (r) and electron-ion RDF g A (r) as follows: 

n\gn(r) = n,(r| I) = n\(r\Uf) = n\ exp[-8Uf(r)] , (2) 
nigAr) = nM I) = nl (rp f) . £ exp| ^ )| + 1 , (3) 

where /ig denotes the chemical potential of a non-interacting electron gas, n\ {n^) is the 
number density of ions (electrons), and (3 = (k^T) -1 the inverse temperature. The electron- 
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density distribution n®(r\U) is determined by solving the wave equation for an electron under 
the external potential U (r) 

-(h 2 /2m)V 2 + U(r)] ^(r) = e^r) . (4) 

In the similar way to the case of classical binary mixtures, the effective external potentials 
Uf G (r) given by Eq. (|T|) are written as 

Uf{r) = v a (r) - r a (r)/P - B a (r)/(3 , (5) 

r a (r) = W C a (\v - r'\)n l [g a (r) - l]dv' , (6) 
i J 

in terms of the direct correlation functions (DCF's) Cij(r) and the bridge functions Bn(r). 
Here, the DCF's Cij(r) in the ion-electron mixture are defined within the framework of the 



DF theory by 

,£ 2 .Fint[n/,7i e 



C^-(|r-r'j) = -/3- 



(7) 



5ni(r)5rij(r') 

where the suffix denotes the functional derivative at the uniform densities 0. Actually 
the explicit expression for the DCF's are given by the Fourier transform in the matrix form 

vifc(k)Vjs = (x^r 1 - (xq)- 1 (8) 

in terms of the density response functions, xq =11 Xy(^) II an d Xq =11 X° l (Q)^ij II; °f the 
interacting and noninteracting systems, respectively, with M =|| n l Sij \\ 0. Note here that 
the density-density response functions Xn(Q) concerning ion reduce to the structure factors 
Sq(Q) an d X 0I {Q) = 1) since the ions behave as classical particles 0. ^From this definition 
of the DCF's the Ornstein-Zernike relations are derived for the ion-electron mixture: 

9n(r) = C„(r) + r„(r) , (9) 
g*(r) = BC a (r)+Br A (r) . (10) 

Here, B denotes an operator defined by 

F Q [By(r)} = (x Q yF Q if(r)} = (x° Q rJ exp[iQ • r]/(r)dr , (11) 



for an arbitrary real number 7, and represents a quantum-effect of the electron through the 
density response function Xq of the non-interacting electron gas. 

A set of integral equations (0) and (|3]) are exact but formal expressions, as well as all other 
equations in the above. However, the ionic charge Z\ and the electron-ion interaction v e i(r) 
must be given beforehand, when we apply these formula to a liquid metal as an ion-electron 
mixture. In order to determine these quantities from first principle, a liquid metal must 
be treated more fundamentally as a mixture of nuclei and electrons (the nucleus-electron 
model), where all interactions between particles are known as pure Coulombic. In this 
model, input data in dealing with a liquid metal is only the atomic number Z\ to specify 
the material. For this purpose, let us consider a liquid metal as a mixture of Nj nuclei 
and Zp^Ni electrons, and solve the problem to determine the electron-density distribution 
around a nucleus fixed at the origin in this mixture. Since a fixed nucleus causes an external 
potential U(r) = —Z^e 2 /r for this mixture to produce an inhomogeneous system, the DF 
theory can be applied to this problem. It should be noticed that DF theory contains some 
arbitrariness in the choice of a reference system to describe the system [[|. We can get a 
simple description of the nucleus-electron mixture if the reference system is chosen to be a 
mixture consisting of JVj — 1 noninteracting ions and Zi(N\ — l) + Z^ noninteracting electrons: 
here, each ion is assumed to have Z-q bound-electrons with a charge distribution pb(r) around 
it and an ionic charge Z\ = Z^ — Z-q. With use of this reference system, the DF theory can 
provide an effective external potential v^(r) for electrons around the fixed nucleus. Then, 
the electron-density distribution n e (r|N) around the fixed nucleus is obtained by solving the 
wave equation for v^(r) m the sum °f the bound- and free-electron parts: 

n e (r|N) = n° b (r|<£) = n c b (r|N) + <(r|N) . (12) 

Hence, the bound-electron distribution n° (r|N) thus determined constitutes the definition of 
the "ion" in the electron-ion model. Furthermore, this bound-electron distribution ng(r|N) 
should be taken identical with the electron distribution Pb[r) of an ion in the reference 
system, since the ion formed around the central nucleus is necessary to be the same structure 
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as any ion in the system. Thus, we obtain a self-consistent condition to determine the 
distribution Pb( r ) in the premise: 

p b (r) = n e b (r|N) ee n e b (r) (13) 

with the bound-electron number = J p^{r)dr. On the other hand, the free-electron part 
n e (r\N) in Eq. fll2|) is taken as the electron-ion RDF n e Q g cl {r) of the electron-ion mixture 
with the free-electron density Uq = Zin l , and the nucleus-nucleus RDF becomes the ion-ion 
RDF g n (r). 

With use of this reference system, we can obtain a tractable expression of f^ 7 ") f° r the 
wave equation to determine n e (r\N) by introducing some approximations to the exchange- 
correlation term involved in it [HI: 



+ J Vee (\r ~ r'\)n h c (r')dr> + pxcKV) + n e ) - MxC K) . (15) 



= Ur) ~ i E / " AWMr') - l]dv> , (14) 

where /ixc(^) is the exchange-correlation potential in the local-density approximation 
(LDA). Note that this expression is equal to Eq. (|5|) without the electron-ion bridge functions 
B e i(r) except that the bare electron- ion interaction is explicitly given by 

_ , v Z A e 2 

v e i{r) = 

r 

In this way, the treatment of a liquid metal as a nucleus-ion mixture is shown to provide 
the ion-electron model, where the bare electron- ion interaction v e i(r) and the ionic structure 
Pb( r ) can be determined in a self-consistent manner. 

With the help of the result from the nucleus-electron model, we can derive a closed set of 
integral equations for the ion-electron mixture, if we introduce the following approximations: 
(A) the electron-ion bridge function in Eq. (|21| ) is neglected: -B e i( r ) — 0, (B) the electron- 
electron DCF C ee (r) is approximated as 

C ee (Q) = ~Pv cc (Q) [l - G jcll (Q)] • (16) 

using the local-field correction (LFC) G^ cll (Q) of the jellium model for an electron gas, (C) 
the bare ion-ion potential i>n(r) is taken as pure Coulombic, i.e., fn(r) = (Zie) 2 /r, and 



(D) the bare electron- ion potential is given by v e i(r) = v e i(r) of Eq. fll5|) . We have called 
this set of equations the quantal hyper-netted chain equations because of the approximation 
B eI (r) ~ in Eq. (|). 



III. MOLECULAR DYNAMICS SIMULATION BASED ON QUANTAL 
HYPER-NETTED CHAIN THEORY 

It is important to realize that the electron-ion model leads to the neutral-fluid model, 
where the ionic behavior of a liquid metal is taken to be the same as a neutral one-component 
fluid interacting via a binary effective interaction in treating the ion-ion RDF. This neutral- 
fluid model is derived from the electron-ion model, when an effective ion-ion potential is 
defined in such a way that the RDF of a one-component fluid should become identical with 
Qu{r) of the electron-ion mixture: 

g(r) = exp[-[3v cS {r) + r(r) + B{r)} = g u {r) , (17) 

with use of the Ornstein-Zernike relation for a neutral one-component fluid 

g {r) - 1 = C(r) + r(r) . (18) 

In the above, the DCF for the one-component fluid is defined by n l C(Q) = 1 — SfiW) 1 
and r(r) = JC(\r - r'\)n l [g(r') - l}dr'. 

Thus, we can write the explicit expression for the effective ion-ion potential of a liquid 
metal in the neutral-fluid model as 



\C el {Q)\ 2 nlx° Q 
l-n e C ee (Q)x° Q 

by taking the bridge function B(r) to be Bu(r) of the electron- ion mixture. Equation ( [3~9[ ) 
can be interpreted within the scope of the standard pseudopotential theory by regarding 
C e i(r) as the pseudopotential Wb(Q) = —C c i{Q)/(3. In this way, the ion-electron model 
is reduced exactly to the neutral- fluid model with a binary effective interaction Eq. (|T9"D. 
in which the many-body forces are taken into account in the form of the linear response 
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expression (0), since the nonlinear effect in the electron screening is involved in terms of 
the electron-ion DCF, which plays the role of a nonlinear pseudopotential. 

By noting the above relations ( |TTD - (JT9|) , the exact expressions (0) - (|T0| ) for the electron- 
ion model can be transformed into a set of integral equations: one is the integral equation 
for a one-component fluid with the effective ion- ion potential v e s(r) 

C(r) = exp[-(3v cS (r) + T(r) + B(r)} - 1 - r(r) , (20) 

and the other an equation for the effective ion- ion interaction v e g(r), that is expressed in 
the form of an integral equation for the electron-ion DCF C e i(r): 

BC el (r) = n° e (r\v el - r cl / (3 - B cl /(3)/n c - 1 - Br eI (r) , (21) 

since the effective interaction v e g(r) is given in terms of C e i(r) by Eq. flT9|). In contrast 
with the usual effective potential in the pseudopotential theory, the effective potential (|19D 
depends on the ion configuration represented by the ion-ion RDF g u (r) through the term: 
r e i(r) = Ej/Cddr - r'\)n l [g a (r) - l]dr> in Eq. (|2J. 

Under the assumptions (A) - (D) mentioned before, the QHNC equations fl2"0p and ( PTj ) 
enables us to perform an ab initio molecular-dynamics (MD) simulation which requires only 
the atomic number Za and thermodynamic states as input parameter, in principle. The 
first estimation for v c g(r) can be obtained with use of C e i(r) evaluated by Eq. (|2l|) with 
an initial guess for g u {f)- Next, an integral equation ( p0|) for a one-component fluid can be 
solved by performing the classical MD simulation for this v e e(r) to produce new ion-ion RDF 
g u (r)] this is used again in Eq. ([H]) to determine a new estimation for v e g(r). This process 
will be continued until convergence of the effective ion-ion potential is achieved (we refer to 
this procedure as the QHNC-MD method). However, such a straight-forward repetition of 
the MD simulation to solve the QHNC equations is not practical in the viewpoint of the 
computational cost. Since the dependence of the effective ion-ion potential f e fr(^) on the 
ionic configuration is rather weak in a simple metal as we have shown in [^], we can adopt an 
approximate theory for B(r) in Eq. (|20|) to get an initial v e s(r) for the QHNC-MD method. 
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For this purpose, we take the variational modified HNC (VMHNC) equation proposed by 



Rosenfeld II], in which the bridge function is approximated by BpY(r]rj) of the Percus- 
Yevick equation for hard-spheres of diameter a with the packing fraction rj = 7rn l cr 3 /6. In 
the VMHNC equation, the adjustable parameter rj is determined by the following condition: 

I„;/ [s M- 9 pv(.;,)]«^* + ^^, (22, 

where gpY^', v) is the RDF for the hard-sphere fluid with the Percus-Yevick equation. Thus, 
in a similar way to the QHNC-MD method, the integral equation (|20"D in the VMHNC 
approximation is solved in a coupled manner with Eq. (^l|) producing new effective ion-ion 
interaction [referred to as the QHNC-VM method]. Furthermore, an initial potential v e s(r) 
to this QHNC-VM method can be obtained by approximating g n (r) in Eq. (|2T| ) by the 
step function 9(r — a) with the ion-sphere radius a = (A-Kn\/3)~ 1 ^ . When this final f e ff(r) 
from the QHNC-VM method is used as an input to the QHNC-MD method, the convergent 
result can be obtained by a few repetition of MD simulation. Finally our procedure to solve 
the QHNC equation with the MD simulation (the QHNC-MD method) is summarized as 
the flow chart shown in Figure [TJ For an initial potential v e s(r) given by approximation 
9n( r ) = $( r — °0) the QHNC-VM method in the preparation phase yields a good initial guess 
for the QHNC-MD method. Then the MD simulation is repeatedly performed to achieve 
convergence of v e g(r) in the refinement phase. 

There are two important points to be noticed regarding the MD simulation when applied 
to the QHNC-MD method. One is that the computer simulation provides the RDF g(r) 
only within the half of the side length L of the simulation cell. This causes an unavoidable 
truncation error in calculation of the Fourier transform J~Q[g{r) — 1] to be used in the 
evaluation of Eqs. ([19]) and (|2T|). The second point is that the MD simulation is performed 
inevitably on a truncated potential for a liquid metal whose effective ion-ion potential is 
accompanied by a long-ranged oscillatory tail: the computer simulation may yield different 
RDF's depending on the cutoff radius R c of the potential. Recently we have proposed a 
precise procedure [[l(J to improve these two defects at the same time and to get the RDF in 
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the whole range of distance for the full potential v e g(r). This method can be applied even 
to the small-size simulation result for the truncated potential w c (r): 

. . v eS (r) - v eS (R c ) for r < R c 
u c [r) = < . (23) 

I for r > R c 

As the first step of this procedure, we extract the bridge function from the raw MD RDF 
data. For this purpose, we extend the the raw RDF data of the MD simulation (?md(^), by 
solving an integral equation 

I <?MD(r) for r < R 

9(r) = < , (24) 

I exp[— (3u c (r) + -T( r )] for r > R 

coupled with the Ornstein-Zernike relation, where R is the extrapolating distance (R < L/2). 
At this stage, in order to obtain a reliable bridge function, it is essential to take R as short 



as about 3 to 4 interatomic spacings or simply as R = R c |T0[ and to discard the RDF data 
outside the distance R so as to reduce the statistical noise contained in the raw RDF data. 
Using the extended g(r), the bridge function B MD (r) can be extracted for distances where 
S-MD(r) ^ by 

„ , x ^(r) + \n[g MD (r)\ - T(r) for r < R 

5md(h = < • (25) 

I for r > R 

At the second step to get the RDF for the full potential v e s(r), we solve the integral equation 
([201) for the full potential with use of this -Bmd(t') as an approximation to that of the full 
potential: this can be justified by the fact that the bridge function is not sensitive to the 



long-range part of the potential and becomes very weak for long-range distance [10 . 



IV. APPLICATION TO LIQUID ALKALI METALS 
A. Numerical procedure of QHNC-MD method 



Liquid alkali metals constitute "simple" metals in the sense that the bound electrons 
forming an ion are clearly distinguished from the conduction electrons and the overlap 
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of the core electrons are negligible; the approximations (A)-(D) used in the QHNC-MD 
method become quit good ones for these metals. Therefore, we have applied the QHNC-MD 
method for five liquid alkali metals (Li, Na, K, Rb, and Cs) near the melting point using 
the parameters specified in Table |; the temperature and density have been chosen to be 
compared with experimental data in ]T2|-|iT||. Here, the temperature and density of alkali 



liquids are specified by two dimensionless parameters: the plasma parameter T = (3e 2 /a and 
r s = a/a B in units of the Bohr radius a B , with the average ion-sphere radius a. 

In our application of the QHNC-MD simulation to the alkali liquids, the local-field correc- 
tion of the jellium model in Eq. fll6|) is chosen to be that proposed by Geldart and Vosko [22 



since it has a simple structure and gives a good approximation. In Eq. (|15|), the expression 
given by Gunnarsson and Lundqvist [fj3| is adopted as the LDA for the exchange-correlation 
potential (ixc( n )- 

After the preparation of initial effective potential by the QHNC-VM method, two itera- 
tions in the refinement phase of the QHNC-MD method (Figure |1]) are sufficient to obtain a 
convergent solution for alkali liquid metals; 16,000 particles have been used in the first MD 
run (Run-1) and 4,000 particles for the second MD run (Run-2). The MD simulation has 
been performed over 50,000 time steps for Run-1 and 100,000 time steps for Run-2 with the 
cubic periodic boundary conditions; the temperature of the system is kept constant by the 
isokinetic constraint fllS[| . The equations of motion are integrated by a fifth order differential 
algorithm with the time increment At = 0.0025nQ~ 1 ^ 3 (m I/ 3) with the mass of an ion 



m : : the corresponding real time is shown in Table [p. In each MD simulation of Run-1 and 
Run-2 for iterations, the effective potential is cut at the radius R c located at the node of 
the Friedel oscillation of v e s(r) as shown in Table |I[ All MD simulations have been carried 
out on a vector-parallel processor Monte-4 |2(J at Japan Atomic Energy Research Institute. 
The computational time required for 10,000 steps is about 30 to 50 hours for 16,000 particles 
including the sampling of the RDF. 

The integral equation (El) has been solved by an iterative procedure introduced by Ng 



2lfl to extend the raw MD data g M o(r) to the whole r range: the extending distance R 

13 



in this procedure (|24j) is taken to be R c for the whole cases. The number of grid points 
and step size used in numerical integrations are 1024 points and Ar = 0.025a, respectively. 
Using C(r) obtained by the HNC equation as an initial input function, it takes about 10,000 
iterations to achieve convergence. 

In order to examine both numerical and computational efficiency of the QHNC-MD 
method, we have tested the convergence of the RDF by evaluating following consistency 
measure for g(r): 

A 9i (r) = 9i {r) - g^r) , (26) 
\A gi \^^f\Ag^r 2 dX 2 . (27) 

Here Qi{r) is g u (r) obtained by the z-th MD simulation and <?o( r ) is that obtained by the final 
step of the preparation phase. Figure |] shows the consistency measure (p6|) of the present 
QHNC-MD simulation for the case of liquid Li, yielding the values \Ag\\ = 1.53 x 10 _1 
and \Ag 2 \ = 5.92 x 10~ 3 . It is easily seen that the convergence of g u (r) is very fast; the 
difference of g n (r) between Run-1 and Run-2 is situated almost within the statistical error 
of the sampling of the RDF in the MD simulation: this means accuracy of about 3 to 4 
digits is already achieved in the QHNC-MD calculation of Run-1. In a consistent way to 
the convergence of the RDF, a good convergence of the effective ion- ion potential v e e(r) is 
achieved in Run-1. Similar to the case of liquid Li, a good convergence of g n {r) and v e &(r) 
is also achieved for other liquid metals. It should be noted that the preparation phase of 
the QHNC-MD method largely enhances the convergence of v e s(r). 

Concerning the treatment of the raw MD data for the ion-ion RDF, it should be empha- 



sized that the extension procedure []T0]] to obtain the RDF for the full potential applied to 
the raw MD data is indispensable for the present calculation. This situation is illustrated in 
Fig. [3], where the truncation error in g u (r) due to the use of the cutoff potential in the MD 
simulation is so large that the convergence of the QHNC-MD method will not be attained 
with the raw RDF data. In addition, the extended RDF for the full potential is almost 
identical with the result of VMHNC equation for the same potential (Figure ^) for r > 5a, 
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i.e., after the third peak of RDF. This suggests that the long ranged Friedel oscillation of 
v e &{ r ) typically seen for liquid metals is essential for the detailed structure of the RDF at 
long distances, and it is necessary to include the information of the long-range part of v e s(r) 
into <7n(r) in order to obtain a self-consistent solution of the effective ion- ion potential by 
the QHNC-MD calculation. 

It is concluded that the convergence of the present QHNC-MD simulation is well attained 
from both numerical and computational points of view, helped by a good initial estimation 
from the VMHNC equation in the preparation phase; the extension procedure to obtain the 
full-range RDF for the uncut effective potential is also necessary to get the convergence. 

B. Ion-ion and electron-ion radial distribution functions 

Following the procedure mentioned above, we obtain the effective interatomic interaction, 
the ion-ion and electron-ion RDF's, the charge distribution p(r) of neutral pseudoatom, the 
bound electron distribution n^(r) forming an ion and the bridge functions for liquid alkali 
metals in a self-consistent way from the atomic number as the only input. In the first place, 
we show structure factors, the Fourier transform of the ion-ion RDF's. It is importa nt for 
a detail comparison with experiment to use the MD RDF corrected for a full potential and 
extrapolated to whole range of distance in its Fourier transform. 

Figure ^ exhibits the structure factors for liquid Li calculated by the QHNC-MD simula- 
tion in comparison with the experimental result ||12|| . It is clearly seen that the present result 
is in excellent agreement with the experiment, improving the detailed structure of Sn(Q) 
near its second peak compared with the result of the VMHNC equation (the final result 
of the preparation phase). This improvement on Su(Q) by the refinement phase essentially 
relies on the detailed oscillatory behavior of the bridge function extracted from the raw RDF 
of the MD simulation. The discrepancy in the bridge function between the MD simulation 
and VMHNC equation gives no serious effects on the first peak of the RDF. But details of 
the RDF for 2a < r < 5a are rather sensitive to the oscillation of B(r) in a similar way as 
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discussed in |IJ. Therefore the use of the extracted bridge function B MO (r) to determine 
the corrected RDF is important in order to guarantee the correct behavior of the structure 
factor by the Fourier transform. The structure factors calculated for Na and K are shown 
in Fig. 5 in comparison with the neutron and X-ray experiments |]T3| , [T4|] . The QHNC-MD 
structure factor of Na is in excellent agreement with the X-ray result (full circles) showing 
a small deviation from the neutron data (circles) at the first peak, while the curves of the 
QHNC-MD structure factor for K lies between the neutron (circles) and X-ray (full circles) 
experimental results. Also, the structure factors from the QHNC-MD method for Rb and Cs 
are compared with the neutron and X-ray experiments [j3~5| in Fig. 6. The first peak of Rb 
structure factor observed by neutron experiment (full circles) [[15|] is shifted a little to large 
Q side compared with that observed the X-ray experiment (circles) [[17]]; the QHNC-MD 
result has the same first-peak position to the X-ray data with a little different height and 
shows overall agreement with the neutron observation. On the other hand, the QHNC-MD 
structure factor of Cs becomes higher in the first-peak than the experiment |13| and shows a 
small deviation in the phase of oscillation of structure factor for the large Q region compared 
with experiment. In our treatment, we solve the Schrodinger equation to obtain the bound- 
electron distribution and electron-ion RDF; the relativistic effect is not taken into account. 
In the case of Cs, the relativistic effect may contribute to the calculation of the structure 
factor, since its atomic number = 55 is rather large; there is a possibility that this effect 
may be ascribed to this small discrepancy between the calculated and experimental results 
in Cs structure factor. In the same way as shown in the Li structure factor, the QHNC-VM 
structure factors deviate from the QHNC-MD results near the their second peak for Na, 
K, Rb and Cs as seen Figs. 5 and 6. This second-peak difference between the QHNC-MD 
and QHNC-VM structure factors brings about a distinct difference in the RDF's obtained 
from their Fourier transform, as can be seen from an example of Na shown in Fig. 7. Thus, 
the RDF from the VMHNC equation is found not so exact as to compare with details of 
experimental results for alkali liquids; the QHNC-MD method is shown to produce reliable 
results for all alkali liquids. 
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Next, we proceed to discuss the electronic structure around ion. The electron-ion RDF's 
from the QHNC-MD method are shown in Fig. 8 for the case of Rb at temperature 313 K 
together with the ion-ion RDF. The electron- ion structure factor S c i(Q) is represented in 
the form 0: 

S e i(Q) = ^&ii(Q) , (28) 
in terms of the ion-ion structure factor and the charge distribution of the pseudoatom: 

P{Q > - 1 - n§CUQ)4 ' (2J) 
With the help of the above equations, the pseudopotential method using the Ashcroft model 



potential can evaluate the electron-ion RDF by inserting it in (p9|) instead of —C e i(r)/f3; its 
result using the empty core radius 1.27A has no inner-core structure near the origin as is 
expected. In the QHNC-MD method, we need not introduce any pseudization in treating 
the core region. Therefore, the electron-ion RDF g eI (r) exhibits the inner-core structure 
similar to a free-atom as is shown in Fig. 8 in contrast with result of the pseudopotential 
method (full circles). It is interesting to note that the electron-ion RDF has an oscillation 
with an inverse phase to that of the ion-ion RDF, since the electrons are pushed away by an 
ion as a whole in the core region. The charge distribution p(r) of the pseudoatom in a liquid 
Rb has a similar structure to the distribution ps s (r) of the 5s-electrons in the free-atom; 
therefore, the total electron-density distribution n^{r) + p(r) is almost same to that of a free 
atom as was indicated by Ziman. Also, note that the positions of the dips in the electron-ion 
RDF are coincident with those appeared in the 5s— electron charge distribution ps s (r) in a 
free-atom; these dips in the electron-ion RDF reflect the inner structure of an ion in a metal. 
The electron-ion RDF's of liquid alkali metals (Li, Na, K, Rb and Cs) are shown in Fig. 10, 
where the inner-core structure is omitted near the origin in each curve. The dip in the 
electron-ion RDF becomes shallower and is shifted to the right side indicating the growth 
of the core size, as the atomic number increases from Na to Cs: the curve of Li shows an 
exception to this tendency. 
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When the electron-ion RDF is determined by solving the wave equation for the self- 
consistent potential v^f (r), we obtain from Eq. ( pTj ) the electron- ion DCF C e i(r), which 
yields an effective ion-ion potential f e ff(r) of Eq. (|T9|). Corresponding to each curves of 
g e i{r) in Fig. 10, the effective ion-ion potential v e g(r) is determined for each element as 
shown in Fig. 11. It should be emphasized that there is no units for scaling lengths and 
energies to effective ion-ion potentials determined by the QHNC-MD method in contrast with 
those obtained by using the Ashcroft model potential f24| . Similarly, no scaling features are 



found in the effective ion-ion potential obtained by the Dagens-Rasolt-Taylor method [25 



which can be thought of as an approximation to the QHNC formulation in the sense that 
the ion-ion RDF is there replaced by the spherical vacancy in Eq. ( |14D [f26[| . Nevertheless, 
it is shown that for liquid alkali metals near the melting point the ion-ion RDF's from the 
QHNC-MD method can be scaled almost into one curve by taking the average ion radius a 
in the units of length as indicated in Fig. 12 except Cs. 

As the summarization of this section, we conclude that the present application of the 
QHNC-MD method to liquid alkali metals provides the ionic structures in excellent agree- 
ment with the experiments within the computational capacity available at present, enabling 
to handle a relatively large system size in the MD simulation in contrast with the usual ah 
initio simulation for the electron-ion mixture. So as to be consistent with the ionic structure 
specified by the ion-ion RDF, the QHNC-MD method is shown to give the electron-ion RDF 
g e i(r), the charge distribution p(r) of a pseudoatom and the density distribution rib(r) of 
the bound electrons forming an ion in a liquid metal at the same time; alternatively this 
electronic structure produces the ion- ion potential f e ff(r) used for the MD simulation. 

V. DISCUSSION 

We have shown that the QHNC-MD formulation provides a very precise description of 
simple liquid metals at any state from the atomic number as the only input: this method 
is proven to yield a first-principles calculation. Our previous calculations with the VMHNC 
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equation, which has been used for input to the QHNC-MD method, generated the structure 
factors with a small but systematic deviation from experiments near the second peak. Now 
our QHNC-MD method is shown to correct this deviation and to yield results in excellent 
agreement with experiments in the whole range of wavenumber Q. Even in alkali metals, 
there is no scaling unit of lengths and energies for the effective ion- ion potentials v e g(r) 
determined from the QHNC-MD method; each potential is different from other and reflects 
a difference in the bound-electron structure of an ion in each metal. The situation is different 
in the case of the effective potential calculated by use of the Ashcroft pseudopotential, which 
leads to give almost single potential curve for all alkali liquids if scaled with the proper unit 



24| . However, it is interesting to note that the ion-ion RDF's can be scaled by the unit of 
a, the average ion-sphere radius, which enables plot the results in almost single curve for 
alkali liquids (Li, Na, K, Rb) near the melting point except the case of Cs. This fact was 



already found in the experimental results for the structure factors 



In our QHNC-MD formulation, the exchange-correlation effect expressed by the LFC 
G(Q) and the LDA ^xc{ n ) are introduced from outside the framework of our formulation: 
those of the jellium model, where the ion distribution is replaced by the positive uniform 
background. The jellium model gives a good description for the electrons in alkali metals; the 
present success of the QHNC-MD method depends on this fact. In addition, the structure 
factors of alkali metals calculated by this method are almost independent of what kind of 
LFC to choose. However, it should be noted that the LFC in the QHNC formulation should 
depend on the ion configuration precisely, since it is defined for the electron-electron DCF 
in the electron-ion mixture. 

To compare the MD structure factors with experiments in detail, it is important to ex- 
trapolate the MD RDF to large distances and to correct errors caused by the cut effective 
ion-ion potential used in the MD simulation. Our extrapolation method |L0] is shown very 
efficient in dealing with the raw MD data for a liquid metal with ion-ion potential accompa- 
nied by a long-range Friedel oscillation; this extrapolation method is indispensable to obtain 
a convergent solution in the QHNC-MD method. In order to obtain so reliable bridge func- 
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tion for the extrapolation of the MD RDF, it is necessary to get a reasonable statistical 
accuracy for evaluation of the RDF; for this purpose, the MD simulation must be performed 
for at least several thousand particles taking about 10 10 to 10 11 samples [[UJ . 



The Car-Parrinello MD (CP-MD) method is based on the same ground to the QHNC- 
MD method: the electron-ion mixture model for liquid metals and the jellium model for 
electrons, which are treated by the DF theory. Also, the bare ion-ion interaction is taken as 
a pure Coulombic in both treatment. However, in the CP-MD method, the bare electron-ion 
interaction is approximated by a pseudopotential, which is introduced from outside of the 
CP-MD formulation; this fact makes a contrast with our QHNC-MD method where it is ob- 
tained self-consistently within the framework of the QHNC formulation and no pseudization 
is necessary in treating the core region. As a consequence, the electron-ion RDF extracted 
from the CP-MD result does not have an inner core structure. It should be emphasized 
that in the QHNC formulation the inner electronic structure around a fixed nucleus, such 
as the electron- ion RDF in the core region and the bound-electron distribution n^(r) of an 
ion in a liquid metal, is determined to be consistent with the outer structure, that is, with 
the surrounding ion and electron configurations. Therefore, the QHNC-MD method can be 
applied to a high density plasma, where a usual pseudopotential can not be constructed due 
to the fact the atomic structure in a highly compressed plasma state is quite different from 
that in the vacuum. While in the CP-MD treatment the electron distribution is determined 
for the multi-ion configuration and the ions are considered to be interacting via many-body 
forces, our approach deals with the one-center problem to determine the electron and ion 
distributions around a fixed ion in a liquid metal. As this connection, it should be kept in 
mind that the electron-ion mixture can be exactly taken as a quasi one-component system 
interacting only via a pairwise interaction in the evaluation of the ion-ion RDF for simple 
liquid metals. 

The advantage of the present method against the CP-MD method based on the usual 
pseudopotential theory can be summarized as follows: (1) The present procedure is capa- 
ble to handle a large system size (~ 10 3 to 10 4 particles) in the MD simulation within 
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the computational resources available at present, helped by the good initial guess with the 
VMHNC approximation for solving the QHNC equation. (2) In the QHNC-MD method, 
the many-body forces and nonlinear effect in the electron screening are taken into account 
automatically in the form of a pairwise interaction in such a way that the nonlinear pseudopo- 
tential is constructed in terms of C e i(r). (3) By setting up an additional integral equation 
for Cee(r), our method can treat the case where the jellium model for the electrons in a 
metal breaks down, that is, where the exchange-correlation effect begins to depend on the 
the ion configuration ||. Therefore, (4) our method is applicable to high density plasmas 
where the ionic structure becomes significantly so different from a free atom due to high 
compression that the usual pseudopotential theory cannot be applied. 

The QHNC-MD method is developed on the base of the ion-electron model, where the 
bound electrons are assumed to be clearly distinguished from the conduction electrons and 
the ions are so rigid and so small that the ion-ion bare interaction is taken as a pure 
Coulombic v u (r) = (Zie) 2 /r. Therefore, our method is applicable only to a simple metallic 
system. In a transition metal, for example, an "ion" cannot be clearly defined since the 
bound electrons are not distinct from the conduction electron, and the overlap of "ions" 
is significant: many-body interactions becomes important. Our method cannot be applied 
to such a case. While the CP-MD method treat electrons in the multi-center problem, our 
method treat them as the single-center problem to determine the density distribution around 
an fixed ion in a liquid metal. Thus, our method cannot describe states of the electrons in a 
multi-center configuration, such as the density of states. Moreover, it should be noticed that 
an orbital of an ion determined by the QHNC-MD simulation may be taken as an average 
orbital in some sense compared with that determined the CP-MD simulation: there, many- 
body forces are changing at every time-step, while the binary ion-ion interaction remains 
constant at every time-step in the QHNC-MD simulation. Already, we have proposed a 
method to treat non-simple metals [[|; the RDF's, the ionic charge and the muffin-tin 
potential are determined by solving a single-center problem with a bare ion-ion interaction 
as an input, while the density of states, the bare ion-ion interaction and the thermodynamic 
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properties are to be obtained as results of the multi-center problem with use of output from 
the single-center problem. 

The CP-MD method can produce the electron-ion RDF and the electron-ion structure 
S e i(Q) is obtained from its Fourier transform. Thus, the charge distribution p(Q) of a 
neutral pseudoatom can be calculated from Eq. Q2"5|) even in the CP-MD simulation, since 
the relation (p8|) between p(Q) and S e i(Q) is exact if a liquid metal can be treated as a 
ion-electron mixture. Also, the electron- ion DCF C e \{Q) is determined from Eq. ( p9l) to give 
an effective ion- ion interaction v e g(r) from Eq. (|19"D ; these quantities must be consistent to 
the ion-ion structure factor Su(Q), if the LFC G(Q) in the electron-electron DCF C ee (r) is 
chosen exact one in the electron-ion mixture. In this connection, there is a point to notice 
regarding the CP-MD method that the exchange-correlation effect of the valence electrons is 
treated by the LDA approximation: the LFC G(Q) does not appear in the CP-MD method. 
This consistency test in the CP-MD method can be exemplified by applying it to liquid 
alkali metals, which are typical examples of the electron-ion model. At the same time, 
these quantities in the CP-MD method can be compared with the results of the QHNC-MD 
method, both methods providing first principles calculations. 
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TABLES 

TABLE I. Parameters used in the present QHNC-MD simulations for liquid alkali metals. 
a = (4-7rnQ/3) -1 / 3 is the ion-sphere radius; L = (3e 2 /a and r s = (47rnQ/3) -1 / 3 are the plasma 
parameter and the electron-sphere radius in units of the Bohr radius a B , respectively. R c \ and R C 2 
are the cutoff length of the effective ion- ion potential v e g(r) in the MD simulation for the first MD 
run and second MD run [see text], respectively. 



element 


T(K) 


r 


r s 


At (fs) 


R cl (a) 


R c2 (a) 


Li 


470 


203.1 


3.308 


0.940 


5.88 


5.88 


Na 


373 


209.1 


4.046 


2.349 


7.06 


7.06 


K 


338 


185.9 


5.024 


3.996 


6.76 


6.76 


Rb 


313 


187.2 


5.388 


6.585 


6.81 


6.80 


Cs 


303 


180.3 


5.781 


8.954 


6.83 


6.84 
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FIGURES 

FIG. 1. Flow chart of the QHNC-MD method. The initial potential v e s(r) is determined by 
approximating gu(r) in Equation (pl|) by the step function 9{r — a) with the ion-sphere radius a. 



FIG. 2. The consistency measure Agi(r) defined by Eq. 26) for liquid Li; dot-dashed curve, 



Ag\ (r); dashed curve, Ag%(r). Solid curve is the final result for the ion-ion RDF gu(r). 

FIG. 3. Comparison of the raw MD RDF for liquid Li with the extended RDF: dashed curve 
is the raw RDF data; solid curve is the extended RDF with the full potential by Eq. ( |20| ) with 
-Bmd(^); ° is the result from the VMHNC equation for the present self-consistent potential. L is 
the side length of the simulation cell; L/2 = 12.794a = 22.40A in the MD simulation for liquid 
Li with 4,000 particles. The discrepancy between the raw data and extended data becomes clear 
for r > 4a, i.e., after the second peak of RDF. On the other hand, the extended RDF becomes 
identical to that of the VMHNC equation for r > 5a, i.e., after the third peak of RDF. 

FIG. 4. The ion- ion static structure factor S\\{Q) for liquid Li: solid curve, the QHNC-MD 
result; dotted curve, result of the final step of the preparation phase (the QHNC-VM result); o, 



experimental result taken from [12|. The deviation of the VMHNC result from the experiment is 



corrected by the QHNC-MD method. 

FIG. 5. The ion-ion static structure factors Sn(Q) for liquid Na (the upper) and K (the lower): 
solid curve, the QHNC-MD result; dotted curve, the QHNC-VM result; o and •, the neutron and 



X-ray experiments taken from [13] and [14|, respectively. 



FIG. 6. The ion-ion static structure factors Su(Q) for liquid Rb (the upper) and Cs (the 
lower): solid curve, the QHNC-MD result; dotted curve, the QHNC-MD result (the final step of 
the preparation phase); o and •, the X-ray and neutron experimental results for Rb taken from 
[17] and [15], respectively. Experiment result for Cs denoted by o is taken from [13]. 
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FIG. 7. The ion-ion radial distribution function <7n(r) for liquid Na: solid curve, the QHNC-MD 



result; dotted curve, the QHNC-VM result; •, experimental result taken from [27|. The VMHNC 
approximation brings about small errors in the amplitude and phase compared with the experi- 
mental RDF. 

FIG. 8. The electron-ion and ion-ion radial distribution functions for liquid Rb: solid curve, 
the QHNC-MD result; dotted curves, the QHNC-VM result; •, The electron-ion RDF derived by 
the use of the Ashcroft model potential. The electron-ion RDF has an inverse phase to the ion-ion 
RDF in their oscillation around unity. 

FIG. 9. The electron charge distribution of a pseudoatom in liquid Rb: solid curve, the 
QHNC-MD result; dashed curve, the charge distribution of the 5s-electrons in a free Rb atom; 
•, The pseudopotential result with use of the Ashcroft model potential. 

FIG. 10. The electron-ion RDF's for liquid alkali metals: inner-core structures near the origin 
are omitted. 

FIG. 11. Effective ion-ion potentials for liquid alkali metals: Scaling properties are not found 
in the alkali potentials. 

FIG. 12. The ion-ion RDF's for liquid alkali metals: The RDF's of Li, Na, K and Rb are written 
almost in a single curve in the units of length a, the ion-sphere radius. The RDF of Cs (the dotted 
curve) shows a small deviation; the dashed curve: Li. The difference among the results for Na, K, 
Rb are indiscernible on this scale. 
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